% In The Name of GOD
% 
% This code has been written for analysis of multilayer lamellar gratings.
% The incident wave is TM polarized.
% 
% In the first section (Inputs) you should determine the structure and
% adjust convergence parameters. Here, an example of defining a structure
% is given.
% 
% Example: 
% % First layer:
% xt(:,1)=[0.25 0.75].';
% epst(:,1)=[1 4 1].';
% dl(1)=d1;
% 
% % Second layer: This layer is a homogenous layer:
% xt(:,2)=[0.25 0.75].'; % for homogenous layer it is not important how to
%                          define "xt(:,2)". The only thing that you should
%                          be careful is its size which must be compatible
%                          with other layers.
% epst(:,2)=[1 1 1].';
% dl(2)=d2;
% 
% % Third layer:
% xt(:,3)=[0.25 0.75].';
% epst(:,3)=[1 3 1].';
% dl(3)=d3;
% 
% 
% Incident plane wave(TM)
%                   \    |
%                    \   |
%                     \  |   Region 1 (eps1)
%                      \ |
%       -     ----      ----      ----      ----
%      |     |4444|    |4444|    |4444|    |4444|
%   d1<      |4444|    |4444|    |4444|    |4444|
%      |     |4444|    |4444|    |4444|    |4444|
%       -     ----      ----      ----      ----   -
%                                                   |
%                             Air                    >d2
%                                                   |
%       -     ----      ----      ----      ----   -
%      |     |3333|    |3333|    |3333|    |3333|
%   d3<      |3333|    |3333|    |3333|    |3333|
%      |     |3333|    |3333|    |3333|    |3333|
%       -     ----      ----      ----      ----
%                    <-------->            <---->
%                 Grating Period(=1)         0.5
%                  
%                   Region 3 (eps3)
% 
% The outputs of this program are diffraction efficiencies which are stored
% in vectors 'DE1' and 'DE3' for reflected and transmitted waves,
% respectively. For example for N = 2 we have:
% DE1 = diffraction efficiency of [2 1 0 -1 -2]'th reflected order
% DE3 = diffraction efficiency of [2 1 0 -1 -2]'th transmitted order
 
%%%%%%%
%INPUTS
%%%%%%%
% clear
% ----------------------
% Convergence parameters
% ----------------------
N=1;        % Truncation order
M=1;        % Number of polynomial terms (M=3 is sufficient).
% global fwhm



index_1= [  1 ] ;

for ind1= 1:length(index_1)
    
    index_2  = index_1(1,ind1);
    
    
    
%     if abc==1
%         break
%     end

% ---------------------------
% Global structure parameters
% ---------------------------


period_1=linspace(  1 ,  1 ,  1 );     


for ww1=1:length(period_1)
      pp=period_1(1,ww1);   % wavelength (um)




angle_1=linspace( 0,  0  ,  1  );     %% 角度  ！！！！！



for ss1=1:length(angle_1)

angle_spr=   0  ;

      
 Audata1=load('TiTi2.txt');          %% 金属材料折射率   
 
%   Audata2=load('Ag (Silver) - Palik (0-2um) Copy 1.txt');          %% 金属材料折射率   
 
 



lamb_1=linspace(   0.4  ,  2   , 300);     %% 入射波长   ！！！微米







 
for llm1=1:length(lamb_1)         %% 可以并行运算
      lam01=lamb_1(1,llm1);   % wavelength (um) 
      
%   NN= length(lamb_1);
%   a=  lam01*1000;
% 
% %   Audata1=load('Au (Gold) - Palik Copy 2.txt');          %% 金属材料折射率 
% %       AuWavelengthData1 = Audata1(llm1,1);
% %       
   
  a=lam01*1000;

%  Audata2=load('Fe.txt');          %% 金属材料折射率


%    AuWavelengthData2 = Audata2(:,1);
%       AuRefractivityData2 = Audata2(:,2)-Audata2(:,3)*1i;   %% 计算的标志金属折射率的虚部要负值 所以修改这样
%        n_index2 = Curvefitting(AuWavelengthData2,AuRefractivityData2,a,2);    %% 金属折射率赋值
       
%         n_metal_index= n_index2;
% 
%  Audata1=load('FeFe.txt');   

% n_index1 = Audata1(llm1)   ;   %% 计算的标志金属折射率的虚部要负值 所以修改这样



n_index1 = Audata1(llm1,2)+Audata1(llm1,3)*1i;   %% 计算的标志金属折射率的虚部要负值 所以修改这样
      


% nn(1,llm1)=  n_index1 ;



%       n_metal_index1=t5;     %% 上面一层

%       n_metal_index2=    Audata2(llm1,2)-Audata2(llm1,3)*1i;      %% 底下一层
%              n_metal_index2=   n_index1  ;  
        
p=pp;   %% 周期
lambda = lam01/p  ;     % Free space wavelength normalized to grating period
n1 = 1      ;         % Refractive index of incident medium
n3 = 1  ;     % Refractive index of transmission medium
eps1=n1*n1;     %Permittivity of incident medium
eps3=n3*n3;    %Permittivity of transmission medium
alpha=angle_spr*pi/180+1e-6;% Angle of incidence (in radian)
 

%t1= 0.4490 ; %没用到 % ;
 t1=   0.1011 ;
 t2=   0.0510 ;   %%% 单位微米  
 t3=   0.0397    ;
 t4=   0.0165 ;
 t5=    0.1014  ;
 t6=     2 ;


% 0.09698;;0.04453;;;0.0131;;;0.1171;;




%tttttt= t1  ;


d_1=  t1/p  ;
d_2=   t2/p  ;
d_3=  t3/p  ;
d_4=   t4/p  ;
d_5 =  t5/p ;
d_6=  t6/p ;





% d_3=   d3  ;
% d_4=   d4  ;




% ------------------------------
% definition of grating's layers
% ------------------------------
% xt(:,m) are places of transitions in m'th layer and epst(n,m) is the
% permittivity of grating between xt(n-1,m) and xt(n,m). dl is a vector whose
% m'th element is the thickness (Normalized to the grating period) of m's
% layer.
fill=0.5;
w=   p*fill   ;
xx1=  w/p ; 


n_SiO2=  1.46*1.46 ;
n_Al2O3=  1.73*1.73 ;
n_TiO2=  2.41 *2.41  ;



n_Ti = (n_index1)*(n_index1)  ; 


xt=  [ xx1 xx1 xx1 xx1 xx1 xx1]  ;
epst= [   n_SiO2 n_SiO2
    n_Al2O3 n_Al2O3
    n_TiO2 n_TiO2
    n_Ti n_Ti
    n_SiO2 n_SiO2
    n_Ti n_Ti

  
          ].'   ;
      
      
      
      
      
dl=[ d_1 d_2 d_3 d_4 d_5,d_6];





%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
% -------------------
% initial computations
% -------------------
k0=2*pi/lambda;
k1=k0*(eps1^.5);
k3=k0*(eps3^.5);
ux=cos(alpha);
uz=-sin(alpha);
kx = k1 * sin(alpha) - (-N:N) * 2 * pi;
kz1 = sqrt(k1^2-kx.^2);
kz3 = sqrt(k3^2-kx.^2);
for j=1:2*N+1
    if real(kz1(j))>0
       kz1(j)=-kz1(j);
    end
    if imag(kz1(j))<0
       kz1(j)=-kz1(j);
    end
    if real(kz3(j))<0
       kz3(j)=-kz3(j);
    end
    if imag(kz3(j))>0
       kz3(j)=-kz3(j);
    end
end
Kz1=diag(kz1);
Kz3=diag(kz3);
% ------------------
 
Nlayer = length(dl);
for lcounter=1:Nlayer
    d=dl(lcounter);
    if max(epst(:,lcounter))==min(epst(:,lcounter))
        [r11 r12 r21 r22]=RmatrixTM_hom(epst(1,lcounter),d,lambda,alpha,eps1,N);
    else
        [r11 r12 r21 r22]=RmatrixTM_lam(xt(:,lcounter),epst(:,lcounter),d,lambda,alpha,eps1,N,M);
    end
 
%   -------------------------------
%   R matrix Propagation algorithm
%   -------------------------------
    if lcounter==1
         R11=r11;
         R12=r12;
         R21=r21;
         R22=r22;
    else
         Z=inv(r22-R11);
         Rx11=r11-r12*Z*r21;
         Rx12=r12*Z*R12;
         Rx21=-R21*Z*r21;
         Rx22=R22+R21*Z*R12;
         R11=Rx11;R12=Rx12;R21=Rx21;R22=Rx22;
    end
%   -------------------------------
end
 
% ---------------------------------------------------------------
% Finding reflection and transmission coefficients using R matrix
% ---------------------------------------------------------------
GG3=diag((kx.^2+kz3.^2)./kz3);
GG1=diag((kx.^2+kz1.^2)./kz1);
 
 
I = eye(2*N+1);
G11=I-R22*GG1/k0;
G12=-R21*GG3/k0;
G21=-R12*GG1/k0;
G22=I-R11*GG3/k0;
G=[G11 G12;G21 G22];
 
deltai0=zeros(2*N+1,1);
deltai0(N+1)=1;
b=0;
uprim=(k1*cos(alpha)*ux-k1*sin(alpha)*uz);
b(1:2*N+1)=(-ux*I+uprim*R22/k0)*deltai0;
b(2*N+2:4*N+2)=(uprim*R12/k0)*deltai0;
RT=G\b.';
 
Rx=RT(1:2*N+1).';
Tx=RT(2*N+2:4*N+2).';
 
Rz = - kx .* Rx ./ kz1;
Tz = - kx .* Tx ./ kz3;
% ---------------------------------------------------------------
    
% --------------------------------------
% Evaluation of Diffraction efficiencies
% --------------------------------------
DE1=zeros(1,2*N+1);
DE3=zeros(1,2*N+1);
for j=-N:N
    DE1(j+N+1)=-real(kz1(j+N+1)/k1/cos(alpha))*(abs(Rx(j+N+1))^2+abs(Rz(j+N+1))^2);
    DE3(j+N+1)=real(kz3(j+N+1)/k1/cos(alpha))*(abs(Tx(j+N+1))^2+abs(Tz(j+N+1))^2);
end





% RR(ww1,ss1)=  sum(DE1);

% if (   sum(DE1)>1      )          %% 去除NAN
% 
% RR(1,llm1)=1; abc=1;  %% 多一句abc 条件   abc等于1 则改参数错误
% 
% 
% elseif(  isnan( sum(DE1)    ) )
% 
%     RR(1,llm1)=1;abc=1;
% elseif(isinf(sum(DE1) ))
%     
%     RR(1,llm1)=1;abc=1;
% else
%       
% 
% 
% end

RR(1,llm1)=sum(DE1);      %%  计算反射率
AA(1,llm1)= 1- RR(1,llm1); 

fprintf('\n 菜园疯狂跑了 %d 公里',llm1 )


end



end



end





% IndMin=find(diff(sign(diff(RR)))>0)+1;   %获得局部最小值的位置
% 
% [val  peaks]=min(RR);           %% 曲线最小的点
% if size(IndMin,2) >1 
%     abc=1;
%     
% elseif IndMin ~= peaks
%         abc=1;
%     else
% 
% 
% 
% ggg=isempty(IndMin);
% 
% 
% if ggg==1
%     abc=1;
% else
%     
%        
%        yy1= RR(IndMin);
% 
%        yy1= abs(1-yy1);       %% 变成曲线的落差值
% 
%   %%%  这一行注意
% %   yy1=0.4;
%   if yy1<0.6
%       abc=1;
% %       break
%   end
%        
%        
%        
% 
%    tt1=lamb_1(IndMin);
% 
% 
%    zz1(ind1)= tt1;       %% 获取两个峰的位置  峰的位置对应于角度
%    
%    zz2(ind1)=yy1;    %%落差值
% 
%    if  ind1==1
%    fwhm1=  abs(1+RR(IndMin));
%    fwhm2=fwhm1/2;
%    
%   
%       [~, zz3]=min(abs(RR(  :, IndMin:end  )-fwhm2))  ;   %% 右边的峰值
%    
%    zz4=lamb_1(zz3+IndMin-1);
%     
%       [~, zz5]=min(abs(RR(  :, 1:IndMin  )-fwhm2))    ;   %% 左边峰峰值
%    
%    zz6=lamb_1( zz5  );
%    
%    fwhm=abs(zz4-zz6);
% %    fwhm=-fwhm;
% if fwhm >0.45
%     abc=1;
% end
% 
%    end
% 
% 
% end
% 
% 
% end
% 
% if abc==1
%     clear zz1 zz2 fwhm
%     zz1(1)=0;
%      zz1(2)=0;
%      zz1(3)=0;
%  zz1(4)=0;
%  
%   zz2(1)=0;
%   zz2(2)=0;
%  fwhm=-20;
% end

end

% aaa=1;
% % 
% figure,
% [X,Y] = meshgrid(angle_1,period_1*1000 );
% Z=RR;
% mesh(X,Y,Z);
% az = 0;
% el = 90;
% view(az, el);colorbar;
% % figure, plot(lamb_1, RR    )
% 
% savefig('figure 0_1')

figure;
plot(lamb_1, AA);
xlabel('入射光波长 (\mum)');  % 横坐标，单位微米（μm）
ylabel('最佳吸收率');         % 纵坐标

%%a_1= sum(AA)/1000;
